home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
PMAT12
/
PMAT.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1993-01-23
|
24KB
|
895 lines
{
P-Mat - A Turbo Pascal program for Recursive Matrix Algebra
Version: 1.2
Author: Mark Von Tress, Ph.D.
Date: 01/30/93
Copyright(c) Mark Von Tress 1993
DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
FROM USE OF THIS PROGRAM.
}
Unit pmat;
Interface
Type
fp = ^double;
ap = Array[1..1] Of double;
apptr = ^ap;
app = Array[1..1] Of apptr;
vmatrixptr = ^vmatrix;
vmatrix = Object
r, c : integer;
Function m( i, j: integer ): double;
Function mm( i, j: integer ) : fp;
constructor MakeMatrix( vr, vc: integer );
Destructor KillVmatrix;
Procedure Garbage;
Procedure Show( strng: String );
Procedure infomatrix( strng: String );
private
v,vcheck: ^app;
nelems : longint;
onstack : boolean;
Procedure purgevectors;
Procedure allocvectors( rr, cc: integer );
End;
mstackptr = ^mstack;
mstack = Object
constructor InitMatrixStack;
Destructor KillMatrixStack;
Procedure inclevel;
Procedure declevel;
Function ReturnMat: vmatrixptr;
Function DecReturn: vmatrixptr;
Procedure push( Var a: vmatrixptr );
Procedure nrerror( strng: String );
Procedure DumpStack;
private
nummats, level : integer;
next : mstackptr;
mv : vmatrixptr;
Procedure cleanstack( a: vmatrixptr );
Function pop: vmatrixptr;
End;
Var
dispatch : mstackptr;
Function matequals( Var lop: vmatrixptr; rop: vmatrixptr ) : vmatrixptr;
Function ident( n: integer ): vmatrixptr;
Function Mult( A, B: vmatrixptr ):vmatrixptr;
Function emult( a, b: vmatrixptr ):vmatrixptr;
Function scmult( a: double; b: vmatrixptr ): vmatrixptr;
Function multsc( b: vmatrixptr; a: double ): vmatrixptr;
Function divis( a, b: vmatrixptr ):vmatrixptr;
Function scdivis( a: double; b: vmatrixptr ): vmatrixptr;
Function divissc( b: vmatrixptr; a: double ): vmatrixptr;
Function add( a, b: vmatrixptr ): vmatrixptr;
Function scadd( a: double; b: vmatrixptr ): vmatrixptr;
Function addsc( b: vmatrixptr; a: double ): vmatrixptr;
Function neg( A: vmatrixptr ): vmatrixptr;
Function sub( A, B: vmatrixptr ):vmatrixptr;
Function scsub( A: double; B: vmatrixptr ):vmatrixptr;
Function subsc( B: vmatrixptr; A: double ):vmatrixptr;
Function tran( a: vmatrixptr ): vmatrixptr;
Function sweep( A: vmatrixptr; k1, k2: integer ) : vmatrixptr;
Function inv( Amat: vmatrixptr ): vmatrixptr;
Function fill( rr, cc: integer; d: double ):vmatrixptr;
Function submat( a: vmatrixptr; lr, r, lc, c: integer ): vmatrixptr;
Function cv( a, b: vmatrixptr ): vmatrixptr;
Function ch( a, b: vmatrixptr ): vmatrixptr;
Function vecdiag( a: vmatrixptr ): vmatrixptr;
Function reada( fid: String ): vmatrixptr;
Procedure writea( fid: String; a: vmatrixptr; titlestr: String );
Procedure setwid( wid: integer );
Procedure setdec( decimals: integer );
Function getwid: integer;
Function getdec: integer;
Implementation
{///////////////// vmatrix functions ///////////////}
{$R-}
Procedure vmatrix.allocvectors;
Var i: integer;
Begin
if c > 8192 then
dispatch^.nrerror('A matrix cannot have more than 8192 columns');
getmem( v, r * sizeof( apptr ) );
For i := 1 To r Do
getmem( v^[i], c * sizeof( double ) );
End;
Procedure vmatrix.purgevectors;
Var i: integer;
Begin
If v = Nil Then exit;
For i := r Downto 1 Do
freemem( v^[i], c * sizeof( double ) );
freemem( v, r * sizeof( apptr ) );
End;
{$R+}
constructor vmatrix.MakeMatrix;
Begin
r := vr;
c := vc;
onstack := false;
if c > 8192 then
dispatch^.nrerror('A matrix cannot have more than 8192 columns');
nelems := longint( longint( r ) * longint( c ) ) * longint( sizeof( double ) );
If nelems < maxAvail - 256 Then allocvectors( r, c )
Else dispatch^.nrerror( 'cannot allocate matrix' );
vcheck := v;
End; { MakeMatrix }
Destructor vmatrix.KillVmatrix;
Begin
purgevectors;
vcheck := Nil;
End;
Procedure vmatrix.Garbage;
Var errorint : boolean;
Begin
errorint := false;
If vcheck <> v Then errorint := true;
If v = Nil Then errorint := true;
If r < 1 Then errorint := true;
If c < 1 Then errorint := true;
If errorint Then dispatch^.Nrerror( 'garbage' );
End;
{$R-}
Function vmatrix.m;
Begin
{ access a member on the right side of assignment }
If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
dispatch^.nrerror( 'm: index out of range' );
m := v^[i]^[j];
End;
Function vmatrix.mm;
Begin
{ store a matrix element on the left side of assignment }
If ( i < 1 ) Or( j < 1 ) Or( i > r ) Or( j > c ) Then
dispatch^.nrerror( 'mm: index out of range' );
mm := @v^[i]^[j];
End;
{$R+}
Procedure CopySD( Var source, dest: vmatrixptr );
Var
i,j : integer;
Begin
source^.Garbage;
dest^.Garbage;
If source = dest Then exit;
If source^.onstack Then
With dest^ Do Begin
purgevectors;
nelems := source^.nelems;
v := dispatch^.Pop^.v;
r := source^.r;
c := source^.c;
vcheck := v;
source^.v := Nil;
dispose( source, killvmatrix );
exit;
End;
If source^.v = dest^.v Then
With dest^ Do Begin
r := source^.r;
c := source^.c;
allocvectors( r, c );
vcheck := v;
nelems := source^.nelems;
End;
If ( dest^.r <> source^.r ) Or( dest^.c <> source^.c ) Then
With dest^ Do Begin
purgevectors;
r := source^.r;
c := source^.c;
allocvectors( r, c );
vcheck := v;
nelems := source^.nelems;
End;
For i := 1 To dest^.r Do
For j := 1 To dest^.c Do
dest^.mm( i, j )^ := source^.m( i, j );
End;
Function matequals;
Begin
rop^.Garbage;
lop^.Garbage;
copySD( rop, lop );
dispatch^.CleanStack( lop );
matequals := lop;
End;
{////////////////////// stack functions /////////}
constructor mstack.InitMatrixStack;
Begin
nummats := 0;
level := 1;
mv := Nil;
getmem( next, sizeof( mstack ) );
With next^ Do Begin
level := 0;
next := Nil;
mv := Nil;
nummats := 0;
End;
End;
Destructor mstack.KillMatrixStack;
Begin
freemem( next, sizeof( mstack ) );
End;
Procedure mstack.nrerror;
Begin
writeln( 'Fatal Error: ', strng );
writeln( 'Exiting to system' );
halt;
End;
Procedure mstack.inclevel;
Begin
level := level + 1;
End;
Procedure mstack.declevel;
Begin
level := level - 1;
If level <= 0 Then nrerror( 'declevel: decremented too often' );
End;
Function mstack.ReturnMat;
Begin
ReturnMat := next^.mv;
End;
Function mstack.DecReturn;
Begin
declevel;
DecReturn := next^.mv;
End;
Procedure mstack.push;
Var
t : mstackptr;
Begin
a^.Garbage;
getmem( t, sizeof( mstack ) );
t^.level := dispatch^.level;
t^.next := dispatch^.next;
t^.mv := a;
a^.onstack := true;
dispatch^.next := t;
dispatch^.nummats := dispatch^.nummats + 1;
t^.nummats := dispatch^.nummats;
End;
Function mstack.pop;
Var t : mstackptr;
vv : vmatrixptr;
Begin
If dispatch^.level < 0 Then
nrerror( 'pop: missmatched inclevel-declevel statments, level < 0' );
If ( dispatch^.next^.next = Nil ) Or( Nil = dispatch^.next^.mv ) Then
nrerror( 'pop: popping off the tail' );
t := dispatch^.next;
dispatch^.next := dispatch^.next^.next;
dispatch^.nummats := dispatch^.nummats - 1;
t^.mv^.Garbage;
vv := t^.mv;
freemem( t, sizeof( mstack ) );
t := Nil;
vv^.onstack := false;
pop := vv;
End;
Procedure mstack.cleanstack;
Var vv : vmatrixptr;
Begin
While ( dispatch^.next^.next <> Nil ) And
( dispatch^.next^.level >= dispatch^.level ) Do Begin
vv := pop;
If vv <> a Then
dispose( vv, KillVmatrix );
End;
End;
Procedure mstack.dumpstack;
Var mm : mstackptr;
Begin
mm := dispatch;
writeln( '------------ Stack Dump --------------' );
While mm <> Nil Do Begin
With mm^ Do Begin
writeln( '----------------------------------------------' );
writeln( 'nummats, level, next, vector: ', nummats, ' ', level, ' ',
longint( next ), ' ', longint( mv ) );
If mv <> Nil Then mv^.infomatrix( 'stack matrix information' )
Else
writeln;
End;
mm := mm^.next;
End;
End;
{////////////////// matrix functions ///////////////////}
Procedure vmatrix.Show;
Var i,j,w,d: integer;
Begin
w := getwid;
d := getdec;
writeln;
writeln( strng );
writeln;
For i := 1 To r Do Begin
For j := 1 To c Do
write( m( i, j ): w: d, ' ' );
writeln;
End;
writeln;
End; { show }
Procedure vmatrix.infomatrix;
Begin
writeln;
writeln( strng );
writeln;
writeln( 'r,c,size : ', r, ' ', c, ' ', nelems, ' bytes' );
writeln( 'v, vcheck: ', longint( v ), ' ', longint( vcheck ) );
End; { show }
Function mult;
Var i,j,k: integer;
C: vmatrixptr;
s: extended;
Begin
a^.Garbage;
b^.Garbage;
If A^.c <> B^.r Then
dispatch^.nrerror( 'mult: matrices do not conform' )
Else Begin
new( C, MakeMatrix( a^.r, b^.c ) );
For i := 1 To A^.r Do
For j := 1 To B^.c Do Begin
s := 0;
For k := 1 To B^.r Do
s := s + A^.m( i, k ) * B^.m( k, j );
C^.mm( i, j )^ := s;
End;
End;
Dispatch^.Push( C );
Mult := Dispatch^.ReturnMat;
End; {mult}
Function Emult{( A, B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
a^.Garbage;
b^.Garbage;
If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
dispatch^.Nrerror( 'Emult: matrices do not conform' )
Else Begin
new( c, MakeMatrix( A^.r, A^.c ) );
For i := 1 To a^.r Do
For j := 1 To a^.c Do
C^.mm( i, j )^ := A^.m( i, j ) * B^.m( i, j );
End;
Dispatch^.push( c );
emult := Dispatch^.returnmat;
End; {emult}
Function scmult{( a: double; b: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := a * B^.m( i, j );
Dispatch^.push( c );
scmult := Dispatch^.returnmat;
End; {scmult}
Function multsc{( b: vmatrixptr; a: double): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := a * B^.m( i, j );
Dispatch^.push( c );
multsc := Dispatch^.returnmat;
End; {multsc}
Function divis{( A, B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
dd: double;
Begin
a^.Garbage;
b^.Garbage;
If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
dispatch^.Nrerror( 'Emult: matrices do not conform' )
Else Begin
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To a^.r Do
For j := 1 To a^.c Do Begin
dd := b^.m( i, j );
If dd = 0.0 Then
dispatch^.nrerror( 'DIVIS: zero divisor' );
C^.mm( i, j )^ := A^.m( i, j ) / dd;
End;
End;
Dispatch^.push( c );
divis := Dispatch^.returnmat;
End; {divis}
Function scdivis{( A : double; B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
dd: double;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do Begin
dd := b^.m( i, j );
If dd = 0.0 Then
dispatch^.nrerror( 'SCDIVIS: zero divisor' );
C^.mm( i, j )^ := A / dd;
End;
Dispatch^.push( c );
scdivis := Dispatch^.returnmat;
End; {divis}
Function divissc{( A : double; B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
If a = 0.0 Then
dispatch^.nrerror( 'DIVISSC: zero divisor' );
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do Begin
C^.mm( i, j )^ := b^.m( i, j ) / A;
End;
Dispatch^.push( c );
divissc := Dispatch^.returnmat;
End; {divissc}
Function ADD{( A, B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
a^.Garbage;
b^.Garbage;
If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
dispatch^.Nrerror( 'add: matrices do not conform' )
Else Begin
new( c, MakeMatrix( A^.r, A^.c ) );
For i := 1 To a^.r Do
For j := 1 To a^.c Do
C^.mm( i, j )^ := A^.m( i, j ) + B^.m( i, j );
End;
Dispatch^.push( c );
ADD := Dispatch^.returnmat;
End; {add}
Function scADD{( A : double; B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := A + B^.m( i, j );
Dispatch^.push( c );
scADD := Dispatch^.returnmat;
End; {scadd}
Function ADDSC{( B : vmatrixptr; A : double): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := A + B^.m( i, j );
Dispatch^.push( c );
ADDsc := Dispatch^.returnmat;
End; {addsc}
Function neg;
Var i,j: integer;
B: vmatrixptr;
Begin
a^.garbage;
new( b, MakeMatrix( a^.r, a^.c ) );
For i := 1 To a^.r Do
For j := 1 To a^.c Do
B^.mm( i, j )^ := - A^.m( i, j );
dispatch^.Push( b );
neg := dispatch^.returnmat;
End; {negate}
Function sub{( A, B: vmatrixptr): vmatrixptr};
Var i,j: integer;
C: vmatrixptr;
Begin
a^.Garbage;
b^.Garbage;
If ( A^.c <> B^.c ) Or( A^.r <> B^.r ) Then
dispatch^.Nrerror( 'sub: matrices do not conform' )
Else Begin
new( c, MakeMatrix( A^.r, A^.c ) );
For i := 1 To a^.r Do
For j := 1 To a^.c Do
C^.mm( i, j )^ := A^.m( i, j ) - B^.m( i, j );
End;
Dispatch^.push( c );
sub := Dispatch^.returnmat;
End; {sub}
Function scsub{ (a: double; b:vmatrixptr): vmatrixptr;};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := A - B^.m( i, j );
dispatch^.Push( c );
scSub := dispatch^.returnmat;
End; { scsub }
Function subsc{(b: vmatrixptr; a : double): vmatrixptr;};
Var i,j: integer;
C: vmatrixptr;
Begin
b^.Garbage;
new( c, MakeMatrix( b^.r, b^.c ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
C^.mm( i, j )^ := B^.m( i, j ) - A;
dispatch^.Push( c );
Subsc := dispatch^.returnmat;
End; { subsc }
Function tran;
Var i,j: integer;
c : vmatrixptr;
Begin
a^.garbage;
new( c, MakeMatrix( A^.c, A^.r ) );
For i := 1 To A^.r Do
For j := 1 To A^.c Do
c^.mm( j, i )^ := A^.m( i, j );
dispatch^.Push( c );
tran := dispatch^.returnmat;
End; {transpose}
Function sweep;
{
SUBROUTINE TO SWEEP THE NXN MATRIX A ON ITS K1 THRU K2
DIACONAL ELEMENTS (SWP(K)SWP(K)A=A)
INPUT : A Matrix to be swept
K1,K2: elements to sweep; if k1=1 and k2= n then gets inverse
}
Const Eps = 1.0E-8;
Var I,J,k,n,it: integer;
d,ss: extended;
Begin
a^.garbage;
If A^.r <> A^.c Then
dispatch^.nrerror( 'sweep: Amat not square in Sweep: ' );
n := A^.r;
If K2 < K1 Then Begin
k := k1; k1 := k2; k2 := k;
End;
For K := K1 To K2 Do Begin
If Abs( A^.m( K, K ) ) < eps Then Begin
For it := 1 To n Do Begin
a^.mm( it, k )^ := 0;
a^.mm( k, it )^ := 0;
End;
End
Else Begin
D := 1.0 / A^.m( K, K );
A^.mm( K, K )^ := 1.0;
For I := 1 To N Do A^.mm( K, I )^ := D * A^.m( K, i );
For J := 1 To N Do
If J <> K Then A^.mm( J, k )^ :=
- D * A^.m( J, K );
For J := 1 To N Do
If J <> K Then
For I := 1 To N Do Begin
ss := A^.m( J, I );
If I <> K Then A^.mm( J, I )^ :=
ss + (1.0 / D) * A^.m( J, K ) * A^.m( K, I );
End;
End;
End;
dispatch^.push( a );
sweep := dispatch^.returnmat;
End; {sweep}
Function inv;
Var AMat1: vmatrixptr;
Begin
dispatch^.inclevel;
amat^.garbage;
If AMat^.r <> AMat^.c Then
dispatch^.nrerror( 'Matrix not square in Invert ' );
new( amat1, Makematrix( amat^.r, amat^.r ) );
amat1 := matequals( amat1, amat );
amat1^.show( 'amat1' );
AMat1 := matequals( amat1, Sweep( amat1, 1, amat1^.r ) );
dispatch^.push( amat1 );
inv := dispatch^.Decreturn;
End; {Inv}
{//////////////////// patterned matrices ///////////////////}
Function ident;
Var i,j: integer;
a: vmatrixptr;
Begin
new( a, MakeMatrix( n, n ) );
For i := 1 To n Do
For j := 1 To n Do
If i = j Then a^.mm( i, i )^ := 1.0
Else a^.mm( i, j )^ := 0;
dispatch^.Push( a );
ident := dispatch^.ReturnMat;
End; { ident }
Function fill;
Var i,j: integer;
c : vmatrixptr;
Begin
If ( rr < 1 ) Or( cc < 1 ) Then
dispatch^.nrerror( 'fill: negative rows or columns' );
new( c, makematrix( rr, cc ) ) ;
For i := 1 To rr Do
For j := 1 To cc Do
c^.mm( i, j )^ := d;
dispatch^.push( c );
fill := dispatch^.returnmat;
End;
Function submat;
Var i,j,rc,cc: integer;
b : vmatrixptr;
Begin
a^.Garbage;
If ( lr < 1 ) Or( r > a^.r ) Or( lc < 1 ) Or( c > a^.c ) Then dispatch^.nrerror( 'SUBMAT: index out of range' );
If ( lr > r ) Or( lc > c ) Then
dispatch^.nrerror( 'SUBMAT: indexes out of order, lc>c or lr>r' );
new( b, makematrix( r - lr + 1, c - lc + 1 ) );
rc := b^.r;
cc := b^.c;
For i := 1 To rc Do
For j := 1 To cc Do
b^.mm( i, j )^ := a^.m( lr - 1 + i, lc - 1 + j );
dispatch^.push( b );
submat := dispatch^.ReturnMat;
End;
Function ch;
Var i,j,k: integer;
c : vmatrixptr;
Begin
{ horizontal concatination }
a^.Garbage;
b^.Garbage;
If a^.r <> b^.r Then
dispatch^.nrerror( 'CH: matrices have different number of rows' );
k := a^.c + b^.c;
new( c, makematrix( a^.r, k ) );
For i := 1 To a^.r Do
For j := 1 To k Do
If j <= a^.c Then c^.mm( i, j )^ := a^.m( i, j )
Else c^.mm( i, j )^ := b^.m( i, j - a^.c );
dispatch^.push( c );
ch := dispatch^.ReturnMat;
End; { ch }
Function cv;
Var i,j,k : integer;
c : vmatrixptr;
Begin
{ vertical concatination }
a^.Garbage;
b^.Garbage;
If a^.c <> b^.c Then
dispatch^.nrerror( 'CV: matrices have different number of rows' );
k := a^.r + b^.r;
new( c, MakeMatrix( k, a^.c ) );
For i := 1 To a^.c Do
For j := 1 To k Do
If j <= a^.r Then c^.mm( j, i )^ := a^.m( j, i )
Else c^.mm( j, i )^ := b^.m( j - a^.r, i );
dispatch^.push( c );
cv := dispatch^.ReturnMat;
End; { cv }
Function vecdiag;
Var i,j,r: integer;
b : vmatrixptr;
t : double;
Begin
{ insert a vector into the diagonal of I }
a^.Garbage;
r := 0;
If a^.r = 1 Then r := a^.c;
If a^.c = 1 Then r := a^.r;
If r = 0 Then dispatch^.nrerror( 'vecdiag: a is not a vector' );
new( b, MakeMatrix( r, r ) );
For i := 1 To r Do Begin
If a^.r = 1 Then t := a^.m( 1, i )
Else t := a^.m( i, 1 );
For j := 1 To r Do
If i = j Then b^.mm( i, j )^ := t
Else b^.mm( i, j )^ := 0.0
End;
dispatch^.push( b );
vecdiag := dispatch^.ReturnMat;
End; {vecdiag}
Function FileExists( FileName: String )
: Boolean;
{ Returns True if file exists; otherwise,
it returns False. Closes the file if
it exists. }
Var
f: File;
Begin
{$I-}
Assign( f, FileName );
Reset( f );
Close( f );
{$I+}
FileExists := (IOResult = 0) And
(FileName <> '');
End; { FileExists }
Function ReadA;
Var i,j : integer;
b : vmatrixptr;
f : text;
titlestr: String;
Begin
If Not FileExists( fid ) Then Begin
writeln( 'trying to open : ', fid );
dispatch^.nrerror( 'reada: file does not exist' );
End;
assign( f, fid );
reset( f );
writeln( 'reading: ', fid );
readln( f, titlestr );
writeln( 'title:' );
writeln( titlestr );
read( f, i, j );
new( b, makematrix( i, j ) );
For i := 1 To b^.r Do
For j := 1 To b^.c Do
read( f, b^.mm( i, j )^ );
close( f );
dispatch^.push( b );
reada := dispatch^.returnmat;
End; { reada }
Procedure writea;
Var i,j,w,d : integer;
f : text;
Begin
assign( f, fid );
rewrite( f );
w := getwid;
d := getdec;
writeln( f, titlestr );
writeln( f, a^.r, ' ', a^.c );
For i := 1 To a^.r Do Begin
For j := 1 To a^.c Do
write( f, a^.m( i, j ): w: d, ' ' );
writeln( f );
End;
close( f );
End; { reada }
Var
printdec, printwid : integer;
Function getwid;
Begin
getwid := printwid;
End;
Function getdec;
Begin
getdec := printdec;
End;
Procedure setwid;
Begin
If wid > 0 Then printwid := wid;
End;
Procedure setdec;
Begin
If ( decimals > 0 ) And( decimals <= printwid ) Then
printdec := decimals;
End;
Begin
setwid( 6 );
setdec( 2 );
new( dispatch, InitMatrixStack );
End.